#!/usr/bin/python
"""
This example demonstrates the FFT of a simple sine wave and displays its
bilateral spectrum.  Since the frequency of the sine wave is folded by
whole number freqStep, the bilateral spectrum will display two non-zero point.
 
Note:
 
This example is coded original in Matlab from Roger Jang's
Audio Signal Processing page.  I translated it into Python with matplotlib.
 
See Also:
 
- "Discrete Fourier Transform" by Roger Jang
    <http://140.114.76.148/jang/books/audioSignalProcessing/ftDiscrete.asp>
"""
__author__ = "Jiang Yu-Kuan, yukuan.jiang(at)gmail.com"
__date__ = "December 2006"
__revision__ = "1.1"
 
import math
import numpy
import pylab as pl
import generation

if __name__ == '__main__':
 
#def fftshift(X):
#    """Shift zero-frequency component to center of spectrum.
 
#    Y = fftshift(X) rearranges the outputs of fft
#    by moving the zero-frequency component to the center of the array.
#    """
#    Y = X.copy()
#    Y[:N/2], Y[N/2:] = X[N/2:], X[:N/2]
#    return Y
 
    N = 32              # the number of points
    Fs = 8000.          # the sampling rate
    Ts = 1./Fs          # the sampling period
    freqStep = Fs/N     # resolution of the frequency in frequency domain
    print "freqStep = %f"%freqStep
    f = 3*freqStep     # frequency of the sine wave; folded by integer freqStep
    print f 
    t = numpy.arange(N)*Ts         # x ticks in time domain, t = n*Ts
    y = numpy.cos(2*math.pi*f*t)   # Signal to analyze

    Y = numpy.fft.fft(y)           # Spectrum
    print type(Y)
    print type(Y[0])
#   print Y 
    Y = numpy.fft.fftshift(Y)      # middles the zero-point's axis
     
    pl.figure(figsize=(8,8))
    pl.subplots_adjust(hspace=.4)
     
    # Plot time data
    pl.subplot(3,1,1)
    pl.plot(t, y, '.-')
    pl.grid("on")
    pl.xlabel('Time (seconds)')
    pl.ylabel('Amplitude')
    pl.title('Sinusoidal signals')
    pl.axis('tight')
     
    freq = freqStep * numpy.arange(-N/2, N/2)  # x ticks in frequency domain
     
    # Plot spectral magnitude
    pl.subplot(3,1,2)
    pl.plot(freq, abs(Y), '.-b')
    pl.grid("on")
    pl.xlabel('Frequency')
    pl.ylabel('Magnitude (Linear)')
     
    # Plot phase
    pl.subplot(3,1,3)
    pl.plot(freq, numpy.angle(Y), '.-b')
    pl.grid("on")
    pl.xlabel('Frequency')
    pl.ylabel('Phase (Radian)')
     
    pl.show()